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Higher-order perturbative calculations in Quantum (Field) Theory suffer from the factorial increase 
of the number of individual diagrams. Here I describe an approach which evaluates the total 
contribution numerically for finite temperature from the cumulant expansion of the corresponding 
observable followed by an extrapolation to zero temperature. This method (originally proposed 
by Bogolyubov and Plechko) is applied to the calculation of higher-order terms for the ground- 
state energy of the polaron. Using state-of-the-art multidimensional integration routines two new 
coefficients are obtained corresponding to a four- and five-loop calculation. Several analytical and 
numerical procedures have been implemented which were crucial for obtaining reliable results. 
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I. Introduction 



Highly accurate measurements require precise theoretical calculations which perturbation theory can 
yield if the coupling constant is small. However, in Quantum Field Theory the number of diagrams 
grows factorially with the order of perturbation theory and they become more and more complicated 
as the corresponding loop diagrams involve high-dimensional integrals over complicated (and singular) 
functions. 

The prime example is the anomalous magnetic moment of the electron where new experiments 
[1, 2] need high-order quantum-electrodynamical calculations. In fact, the estimate for the fifth-order 
contribution is the largest source of theoretical uncertainty if one attributes an "error" to it at all [3] . 
In addition, further improvements of the experimental accuracy are foreseen. 

As derived in the textbook [4] the number of diagrams contributing to the vertex function in 
Quantum Electrodynamics (QED) is given by the coefficients of the generating function 



km 

K (z) 



(1) 



[with z = — l/(4a) and Kq(z) the zeroth-order modified Bessel function of second kind] when expanded 
in powers of the fine-structure constant a 

r(a) = l + a + 7a 2 + 72a 3 + 891 a 4 + 12672a 5 + 202770a 6 + ... . (2) 

The contributions up to third order are known analytically [5] and the 891 diagrams in fourth order 
have been evaluated numerically by Kinoshita and coworkers [6]. In view of the ever more precise 
experiments there are ongoing efforts [7] to calculate all 12672 diagrams in 0(a 5 ) numerically and by 
automated routines. This is a huge, heroic effort considering the complexity of individual diagrams, 
the large cancellations among them and the intricacies of infrared and ultraviolet divergencies in the 
integrands. 

Obviously new and more efficient methods would be most welcome for a cross-check as well as 
further progress. However, it is useful first to consider a simpler field theory which is nontrivial but 
free from ultraviolet divergencies. This is supplied by the polaron problem - the field theory of a 
single nonrelativistic electron slowly moving in a polarizable crystal and thereby interacting with an 
infinite number of phonons. Similar as in Quantum Electrodynamics there exists a large number of 
perturbative calculations for the ground-state energy and other properties of the quasiparticle which 
is made up by the electron and its surrounding cloud of virtual phonons. 

In this paper we investigate a method originally proposed by Bogolyubov (Jr.) and Plechko (BP) 
[8] to obtain higher-order terms in the ground-state energy of a polaron without evaluating diagrams. 
As the polaron problem is the prototype of the worldline approach to relativistic Quantum Field 
Theory [9, 10, 11, 12] we believe that a similar method also holds promise for high-order perturbative 
calculations in particle physics, in particular QED. 

Preliminary results have already been presented in Ref. [13]. Here I give a detailed account of 
the analytical and numerical methods which are required so that the BP method works. The paper is 
organized as follows. In Sees. II. and III. we recall the basics of the polaron model and the BP method. 
Section IV. gives an account of the necessary steps to obtain reliable numerical results. These are 
presented and discussed in Sec. V.. The last section contains our conclusion and the outlook for 
further work whereas more technical details are collected in three appendices. 
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II. The polaron problem - a nonrelativistic field theory 

A model Hamiltonian describing the dressing of the bare electron by a cloud of phonons has been 
given by H. Frohlich 

H = \v 2 + J d 3 kala k + i(2V27Ta) 1/2 J ^ ^ [4 e*"* " H.c. } , [a k ,a£/] = ^(k-k') (3) 

where a is the dimensionless electron-phonon coupling constant. Due to its interaction with the 
medium the energy of the quasiparticle is changed and it acquires an effective mass 

The aim is to calculate the power series expansion for the ground-state energy of a non-moving polaron 

E (a) =:Y / e n a n (5) 

n=l 

as function of a [14]. The lowest-order coefficients are well known 

ei = -1 (6) 
e 2 = - In (\ + ^y/2j = -0.015919622 from Ref. [15] , (7) 
e 3 = -0.000806070 from refs. [16, 17] , (8) 

but there has been no progress towards higher-order terms. 

In the path-integral approach [18] the (infinite) phonon degrees of freedom may be integrated out 
exactly which leads to an effective, two-time action 

LJ Jo 2 2V2 Jo sinh(/3/2) |x(t)-x(t')| 

Here f3 is the Euclidean time or inverse temperature. Some simplifications are possible: first, the 
symmetry between the two times t, t' allows us to restrict the integration range of the latter to 
< t' < t together with doubling the strength of the interaction. Second, as we are only interested in 
the ground-state energy i?o of the polaron which can be obtained by the large-/? limit of the partition 
function 

Z := [d 3 x[ V 3 xe- s W ^ const e^ E °, (10) 

J Jx(0)=x(/3)=x 

we may replace 

cosh \p/2 - (t - t')] exp(-<r) + exp[-(/?-<7)] /^oc 

— — = — i-ex P (-/?) > eM ~ a) ' 

where a = t — t' is the relative time [20]. Thus, in the following, we will use 

S[x] = f dt\ x 2 - " f dt f* dt'-^t^- = : So + S, (12) 
L J Jo 2 y/2 Jo Jo |x(t)-x(tf)| v ; 

as a full polaron action. 
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Useful order-of-magnitude estimates for higher-order energy coefficients can be obtained in various 
approximate treatements of the polaron problem. Most prominent and successful among these is 
Feynman's approach [18] in which a quadratic trial action 

S t = f dt-± 2 + f dt f dt' f(t-t') [x(t) -x(t')] 2 (13) 
Jo 2 jo jo 

is used as variational approximation for the full action (12). Feynman chose an exponential form of the 
retardation function with two variational parameters which are determined by minimizing Jensen's 
inequality. The corresponding energy coefficients can be calculated analytically to high order [17] as 
sketched in Appendix Appendix A:. The result is 

1 16 56 

ef = -1 , ef = = -1.234568 x 1(T 2 , ef = - \fl = -0.634366 x 1(T 3 , 

1 ' 2 81 3 729 6561 

= 3200/15-633236 78m ^ = 

1594323 531441 ' y J 



e 4 
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1673496632 - 6044800^/10- 70304^/13 793600 r— 1476371144 /- 

1 V70 V7 

129140163 43046721 301327047 

= -0.395686 x 10" 5 . (15) 

However, one can do better by allowing the variational principle to determine the best retardation 
function itself. Then one gets [19, 21] 

e\ cst = -1 , e\ est = - (j^ - 7^) = "1-2597803 x 10~ 2 . (16) 

Note that e2 est is only slightly better than ef despite the fact that the retardation function in the un- 
restricted variational approach has quite a different small-time behavior than Feynman's parametriza- 
tion. This is due to the (relative) insensitivity of the polaron energy to small-time dynamics. In this 
respect four-dimensional field theories in the worldline description are quite different, in particular 
realistic, renormalizable ones similar to QED [12]. Appendix Appendix A: also describes how one 
can obtain numerically the higher-order energy coefficients for the best quadratic approximation. We 
have obtained the values 

e best = _ , 64650 x 1Q -3 ^ e best = _Q 46g6 x 1Q -4 ^ gbcst = _ _ 3940 x 1Q -5 ( 17 ) 

which - again - are not very much different from the results using the much simpler Feynman 
parametrization. 



III. The Bogoliubov-Plechko (BP) method 

In order to get the perturbative expansion of Eq{o) we use the cumulant expansion of the partition 
function for large (5 

(18) 

where X n (P) are the cumulants with respect to Si and Zq is the free partition function for a system 
confined in a large volume. 



Z = Zq exp 



E 

n=l 



MP) 
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The cumulants (or semi-invariants) are obtained from the (normalized) moments 

rx(/3)=x 
/x(0)=x 



//•x(/3)=x 
d 3 x / V 3 x e~ So[> 

JxfO~)=x 



(19) 



(here the (3 dependence is suppressed and the normalization constant C is chosen such that tuq = 1) 
via the recursion relation 



An+i — m n +i — E ( u ) ^ k+1 



m n -k 



This is standard and easily proved by differentiating the characteristic function [22] 



*(t) = (e-* 51 ) = £(-)"-m n = exp 



n=0 



. n=l 



with respect to £ in moment and cumulant form 



E 



n! 



■m n+1 = • ^ 



(-*)» 



n! 



An 



+1 



(20) 



(21) 



(22) 



n=0 n=0 

If the moment expansion for $(i) is inserted on the right-hand side one obtains after rearrangement 



n=0 



n=0 



Afc+i rn n ^ k 



for all powers of t which establishes Eq. (20). The first few cumulants are 

Ai = m\ 

A 2 = m 2 — m\ 

A3 = wis — 3 1JI2 mi + 2 mf 

A4 = m4 — 4m3 mi — 3m| + 12 m2 — 6mf 

A5 = ms — 5 m4 m\ — 10 7723 m2 + 20 771,3 m\ + 30 m\ — 60 7712 mf + 24 mf . 
For large (3 we then get the ground-state energy as zero-temperature limit of the free energy 



£° = }™ -o E 



1 







n=l 



n! 



-An(/3) 



(23) 



(24) 
(25) 
(26) 
(27) 
(28) 



(29) 



since the free partition function does not contribute. By construction the nth moment is proportional 
to a n and Eq. (20) (and the examples) show that the cumulants share this properties. Comparing 
with Eq. (5) we see that 

R n+1 l 



hm - \ n {f3) . 



(30) 



The moments m n . We calculate the moments m n by expanding the paths in Fourier components 



x(t) = V^b | + f:^b fcS in(^) , x=:^b 



so that 



so = E b i 

k=0 



(31) 



(32) 
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and the functional integration is over the coefficients b^, k = 0, 1, Writing 

rP r* . ,. „s r d 3 p 1 



we have 



m. 



/ *!•••*»/ <■■■} < ex p - f/ i) - • • • - (*» - '»)] 



"fwji-Jtfg \ exp 



d 3 p„ l 



where 



and 



4(M') = v ' 2 



m=l fc=0 



: fc = 0, 



(O) : = 



/ d 3 6 d 3 6i • • • O(bo, bi . . .) exp [-5 (b , bi . . .)] 



Jdtbotfh... exp[-5 (b ,bi...)] 
is the average with respect to the free action Sq. 

As a Gaussian integral over the b^'s this average can be done easily and one obtains 



Oi 



2 n/2 



m 

I r J3 



exp 



m=l 



n 

m=l 



2tT 2 P 2 



exp 



od / n 



fc=0 \m=l 



If we now write the mth Coulomb propagator as 



1 1 

Pi " 2 



d-u m exp 



then all momentum integrations can be performed and give the result 



m r 



(-r 



a 



n n 



f/3 rtm roo \ 

(47r) n/2 ll A y dt mJ </ du m je W 



n 



m=l 



x [det n A (t±, . . . , t n , t l5 . . . , t n ; ui, . . . , u n )] ^ 
Here (A) is the n x n matrix made up by the elements 



(33) 



(34) 



(35) 



(36) 



(37) 



(38) 



(39) 



(A)ij = 2 £ k(U, t'i)£ k (tj, t'j) + Ui Sij = : Oij + Ui Sij . (40) 

fc=0 

It is essential that the infinite sum over the modes k can be performed analytically Using Eq. 1.443.3 
in Ref. [23] 

OO ; 2 2 

■c-^ COS KX 7T 7TX X , , . 

E ""Or" = —- — + — > < x < 2^ (41) 



fc=i 



fc 2 



6 
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we indeed have 



^ 1 . kirx . kiry 



1 00 1 
1 



kir(x — y) kir(x + y) 

cos — ^— — cos - 



2/3 



xy 

min(x,y) - — 



< x,y < 



and, therefore, 



aij = min(tj, ij) — min(ij, i^-) — min(t^, tj) + min(^, ^) . 
Using min(x, y) = [x + y — \x — y\] /2 this may also be written as 

1 



Note that a^- = dji and that 



*i - t \ = : Gi > 



(42) 
(43) 

(44) 
(45) 



since t i > t\ . This is a special case of the more general fact that (A) is a positive definite matrix 
(otherwise the momentum integral would not converge) [24] . Well-known theorems of matrix analysis 
[25, 26] then guarantee that the principal minors of all orders are non- negative and the diagonal 
elements are just the ones of lowest order. 

Introducing total and relative times 



Oi : 



ti , £j 



U + 1' 



(46) 



we have 



a n ™ / CP rp-a m /-z poo \ 

II [J Q dX m J Q du m j exp 



(-Y 



P-<Tm/2 



m=l 

x [det n A (cti, . . . , a n , £1, . . . , £ n ; m, . . . , u n )}~ 3/2 



(47) 



Due to time-translational invariance, the nondiagonal matrix elements, say 012, only depend on three 
variables which we denote by 



S : = £1 - £ 2 , r : = ^ (<7i - 02) , s : = ^ (ai + a 2 ) > 



(48) 



Then one has 



Ol2 



|5 + s| + |5 — s\ — \S + r\ — \S — r\ 



s — \r\ for \S\ < \r\ 
s - \S\ for \r\ < \S\ < s (49) 
for \S\>s. 



Figure 1 shows that 012 is indeed a nonanalytic function of the times as expected from the absolute 
values in Eq. (44). Note that it is even in S,r,s. If we would split up the integration region into 
subregions where the time differences have definite sign we would get rid of that complication at the 
price of considering many different contributions. This is exactly what happens in the diagrammatic 
approach and is the source of the proliferation of diagrams in high-order perturbation theory. 
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Figure 1: (Color online) The nondiagonal matrix element a\2 of the matrix (A) as a function of the 
time variables defined in Eq. (48). 



IV. How to make the BP approach numerically feasible 

Equation (30) together with Eqs. (20) and (39) specify how to calculate the nth-order coefficient e n 
for the perturbative expansion of the polaron ground-state energy. Taken at face value one needs to 
evaluate a 3n-dimensional integral at large (asymptotic) values of the inverse temperature (3. While 
this seems doable in principle, it is clear that in practice precise values of e n or the numerical feasibility 
of the whole approach need further improvements and refinements. As these practical questions have 
not been adressed at all in Bogoliubov and Plechko's paper [8] we will describe several steps crucial 
for success. 

A. Additional integrations 

It is obvious that any reduction in the dimensionality of the integral to be evaluated numerically 
will be of great help. As explained above the integrations over the times can only be performed 
by splitting the integration regions in many subregions leading to the time-honored diagrammatic 
approach. However, the dependence on the auxiliary variables Uj is simple and analytic and therefore 
it is possible to perform some of the integrations over them by expanding the n x n determinant 
det n j4 into cofactors [27]. For example, the dependence on u n is simply obtained by expanding with 
respect to the nth row (or column) 

det n A = u n A n + det n A(u n = 0) , (50) 

where A n denotes the determinant of the matrix which is obtained from (A) by removing the nth 
row and the nth column, i. e., it is a special (n — 1) x (n — 1) determinant known as principal minor 
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[28]. Therefore the integration over u n in Eq. (39) can be easily performed: 

V n (l,2,...,n) := / du n d^' 2 A(l,2,...,n) = (51) 
Jo A n *Jdet n A[u n = 0) 

Here we use the short-hand notation i := (tj,iJ,Uj) and the integration over n ra is indicated by 
underlining the nth argument. The dependence on u n _i is obtained similarly: 

A n = u n _i ^ n _i, n + A n (n n _i = 0) , (52) 
det n A(u n = 0) = A n _i(«„ = 0) + det n ,4 (u n _i = u n = 0) . (53) 

Here A n _i >n denotes the determinant (principal minor) of the matrix which is obtained from (A) by 
removing both the (n — l)th and the nth row and column. The subsequent integration over u n _i is 
therefore still an elementary one (u n -\ = u n = is understood in all determinants from now on) 



poo poo 

Pjl,2,..., n- l ,n) := / d«„_i / dn n det" 3/2 A(l, 2, . . . , 

Jo Jo 



f 00 2 1 

= / dn„_i == (54) 

Jo n„_i A n _i in + A n v /u n _ 1 ^4 n _ 1 + det„ ^4 

but depends on the sign of the combination A n A n -\ — A n _i^ n det n A This is fixed since all the 
coefficients in the integrand are principal minors of the positive semidefinite matrix (A) which not 
only are non-negative themselves but also obey the Hadamard-Fischer inequality [Ref. [25], Eq. 7.8.9] 

A n —i A n > A n —\ :n A (55) 

{A = &et n A ). Therefore the integration over u n _i gives [see, e.g., Ref. [29], Eq. 192.11] 



where 



-n n o 1 \ 4 arcsin ^TJJ 

P n (l,2,..., n-l ,n) = — . == -=^= — , (56) 

< XifF := 1-4=^ < 1 (57) 

is non-negative and does not exceed unity as needed for a proper argument of the arcsin function. 

Let us illustrate that for the case n = 2 where all principal minors can be evaluated easily. With 
Eqs. (25) and (39) one then obtains 

A 2 = ^ f dt x dt 2 [ h dt[ t dt' 2 e -(*i+fa-*i-4) [D 2 (l,2)-Di(l)X>i(2)] 

47T Jo Jo Jo 

= a - f d h dt 2 t dt[ t dt> e-<*+»-*-U -J= h ( -7^=) , (58) 
Jo Jo Jo vana22 \van022 / 

where 

f 2 (x) : = arcsi ^) _i, (59) 
x 

Thus the second cumulant (and therefore the second energy coefficient) would vanish without the 
nondiagonal matrix element a±2, i.e., the correlation between the times when the two phonons have 
been emitted(absorbed). 
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B. Extrapolation for (3 — > oo 



A crucial question for the feasibility of the BP approach is how the asymptotic limit [3 = oo is reached. 
From Appendix Appendix B: where the cases n = 1, 2 are treated explicitly we expect 



X n ((3) — /? e „ + d n + 0(e^/^) 
so that from Eq. (30) only a rather slow convergence to the asymptotic value is expected: 



e„ = lim 



e n + 



P 



(60) 



(61) 



This can be greatly improved not by dividing \ n ((3) by (3 but by taking the derivative of A n (/3), i. e., 
considering 

(") n+1 „_ ^09) 



e n (P) : 



lim 



a n n! /3-^oo d(3 
which approaches the asymptotic value exponentially 



e n (P) 



/3^oo 8 

' dp 



[3e n + d n + (e^/VP) = e n + O (e^/^) 



(62) 



(63) 



- at least in the analytical examples given in Appendix Appendix B: for n = 1, 2. 
We therefore will assume that for large enough f3 



e n (P) 



a n -0 

6n + 7P e 



for all values of n in the following. Alternatively, the behavior 

«> + fie-" 



(64) 



(65) 



will be fitted to the numerical data if they are precise enough to determine also the power u n . 

Moreover, evaluating the differentiation with respect to (3 also lowers the dimension of the in- 
tegral which has to be evaluated numerically because the variable (3 enters as upper limits of the 
multidimensional integral (47). Writing the corresponding cumulant as 



: (-)' 



a" 



n 



(4vr)«/2 11 yj > J aj/ 



P-aj/2 



dYij) F n (cri, Si; <72, S 2 ; ...; a n , £ n ) 



(66) 



we find no contribution by differentiating the upper limit of the Oj integration since the range of T,j 
then vanishes. Thus 

dX n 



dp 



a n n ( fP \ n ( rP-Vk/i \ 



Ei=/3-0-;/2 

(67) 



For example, for n = 2 we have 



d\ 2 a 2 ft 3 
dp IT Jq 



\j0-\O2 



-(0-1+0-2) 



0-a 2 /2 
02/2 



cE 2 / 2 



012 



+ (1 2) 



=/3-<ti/2 



(68) 
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C. Symmetrization 

We may exchange simultaneously 



a,-, <-> cr fe ,S fe , j ^ k = 1,... ,n 



(69) 



in the integrand of Eq. (67). There are n! ways of doing that and thus 



5/3 (4vr)«/ 2 



n! 



pcrmut. 



(70) 



Ei=^-<7i/2 



:K ymm (K,S fc }) 



and the domain of integration can be reduced [30]: 



(-a) n fl 3 r ai r n - 1 A -i-r ( fl 3 -^/ 2 \ 

h& L dai L ^ • • I *• ? n (X fc/2 ^ «** 



5 ^ ( 47 V • -« , \->°k/' / E 1= /3- CT ,/2 

(71) 

Again taking n = 2 as simple example we find from Eq. (68) that F| ymm = 2F2 as the integrand is 
already completely symmetric. Hence 



2a 2 ^ r 1 exp(-o-i - a 2 ) 



8X 2 2a 2 [P f°i 
— = / da! da 2 

Of) IT Jo JO 



L 



(T1/2 V ^/<?\a 2 



(1-2) 



Ei=/3-cti/2 



(72) 



where we have used an = Gi. Further evaluation of Eq. (72) is presented in Appendix Appendix B:. 
For n > 2 we have to perform the symmetrization explicitly as the integrations over u n ,u n -i lead to 
a nonsymmetric integrand. 



D. Mapping 

Finally for Monte Carlo integration we need a mapping to bring all integration variables into the 
hypercube [0, 1]. After some experimentation we have chosen 

Ui = (Ti^-lj , i = l,2,...,(n-2) (73) 

and 

(7i = (3s\ , ai = Ui-x s 2 , i = 2,...,n, (74) 

= (p-ai) Si + ±<Ti (75) 

as transformation of the remaining variables. Here all £j, Si, Si £ [0, 1]. Equation (74) removes possible 
square-root singularities which are seen in the examples for n = 1, 2 in Appendix Appendix B: - these 
are integrable analytically but would pose severe problems for numerical integration. More refined 
mappings of the relative times (for example, to include the exponential suppression) have been tried 
but did not result in significant improvements. 
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V. Numerical results 



A. A test: e 3 

We have tested our approach by determining the third order coeffcient 63 which has been calculated 
by Smondyrev [16] with later improvements in accuracy [17]. Table I lists the values of es(($) obtained 
by Monte Carlo integration using the classic VEGAS program [31] with umc = 4.9 x 10 8 function 
calls per iterations. We have used 100 iterations for each f3 value. Thus the total number of function 
calls was 

n ( St = nMC n it = 4.9 x 10 10 . (76) 




3 4 5 6 7 8 9 
P 



Figure 2: (Color online) Monte Carlo results for the derivative of the third cumulant as a function 
of the Euclidean time (inverse temperature) j3. The total number of function calls is denoted by ntot 
and the full (open) circles are the points used (not used) in the fit (see Table II. 



Figure 2 shows that e%{0) monotonically approaches Smondyrev's value with increasing (3. The 
sheer fact that e%{($) converges to a constant value at large (3 is a good signal: individual moments 
m n would behave as (3 n for large values of (5 but the construction of the cumulants takes away all 
these j3 powers except the linear one which contains the information about the ground-state energy. 
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Table I: Third order energy coefficient e^{f5) from the derivative of the the third cumulant as a function 
of the inverse temperature [3. The numerical results were obtained with the Monte Carlo routine 
VEGAS for evaluating the full six-dimensional integral. Numbers in parenthesis are the estimated 
errors in units of the last digit. The last column gives the \ 2 P er degree of freedom (Ndf) monitored 
during the iterations. This should be close to one if the iterations are consistent with each other. 



p 


-e 3 (J3) x 10 3 


X 2 /N DF 


4.0 


0.7474 ( 5) 


0.969 


4.5 


0.7704 ( 7) 


0.876 


5.0 


0.7846 ( 8) 


0.836 


5.5 


0.7934 (10) 


0.837 


6.0 


0.7987 (11) 


0.821 


6.5 


0.8017 (13) 


0.792 


7.0 


0.8033 (15) 


0.768 


7.5 


0.8039 (17) 


0.775 


8.0 


0.8041 (19) 


0.772 



Table II: Extrapolation of e^{(3) to (5 = oo using the data from Table I, the fitting range G 
[Pmin, Pmax] and the fixed power 1/3 = 0.5 in the ansatz (64). The last column gives the x 2 /Ndf of 
the two-parameter fit where Njjp = number of data points - 2. 





Anax 


-e 3 x 10 3 


X 2 /N DF 


4.0 


8.0 


0.8043 (6) 


0.989 


4.5 


8.0 


0.8052 (7) 


0.138 


5.0 


8.0 


0.8056 (8) 


0.048 


5.5 


8.0 


0.8055 (10) 


0.058 
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We have fitted these data with the ansatz (64) which, of course, only holds for asymptotic values 
of p. Therefore the lower limit /3 m i n of the fit range [/?min> /^max] was successively raised until the 
X 2 /Ndf of the fit reached a minimum. This is displayed in Table II. If (3 m \ n is too close to (5 m ax 
the degrees of freedom decrease which should cause the x 2 /Ndf to increase in turn [32]. This fitting 
strategy yielded 

e 3 = -0.8056(8) x 10~ 3 . (77) 
If we allow the more general ansatz (65) we obtain as best fit 

e 3 = -0.8055(6) x 10~ 3 (78) 

and 

Vz = 0.55(3) . (79) 

The above error estimates may be a little bit optimistic since we have taken the VEGAS errors at 
face value. In addition, the power f 3 and the parameter a 3 in the fit function (65) turn out to be 
highly correlated. Nevertheless the exp(— [3)/y/j3 behavior also seems to hold for higher cumulants 
and the extrapolated result is in good agreement with Smondyrev's analytical result (8). The main 
message of this test therefore is that (our implementation of) the BP method is working and able to 
give accurate values for the perturbative expansion of the ground-state energy of a polaron. 



B. A new co efficient: e 4 

When applying the previous approach to the calculation of the first unknown coefficient e 4 an un- 
pleasant outcome is found: as seen in Fig. 3 for a fixed value of (3 = 5 the convergence with the 
number of function calls is very slow. Since the cancellations in the integrand are more severe for the 
large [3 which is needed for determining e 4 only a very rough determination of this coefficient was 
possible in acceptable CPU time. 

Fortunately a solution was found by performing the remaining integrations over Uj ,i = 1,2 by a 
deterministic integration routine. While such an option is not available for the time integrations for 
which the integrand is nondifferentiable (see Fig. 1) it is possible for the integration over the auxiliary 
variables Uj where the dependence is an analytic one [see Eqs. (39, 40)]. 

We have used the powerful tanh-sinh integration procedure [33] which - after a judicious transfor- 
mation of variables - is nothing else than the trapezoidal approximation to the transformed integral 

dx fix) Ri h ~ _ 2 _ ~ £ w k f(^-^ + 

with precalculated abscissae x k and weights w^. Since this quadrature rule seems not to be very well 
known (see, however, Ref. [34]) Appendix Appendix C: gives a short account of its basic features 
together with details of our implementation. Having in mind an application to our multidimensional 
case the convergence rate with the number of function calls 

n t = 2k max + 1 (81) 

is of paramount interest. In the one-dimensional case the error may decrease as fast as exp(— cnt/ In nt) 
[35, 36] depending on the analyticity domain of the transformed function f(x). However, without any 
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Figure 3: (Color online) Convergence of the fourth order coefficient e4(/3 = 5) for a fixed value of 
the inverse temperature as a function of the total numer of function calls n tot . Square (blue) 
points denote the case where the full nine-dimensional integral was evaluated by the Monte Carlo 
routine VEGAS, (red) circles show the result if two of the integrations are done by the deterministic 
tanh-sinh-quadrature rule and the rest stochastically. 



knowledge about that and in a multidimensional application, such an error estimate is of no help and 
we have to test the convergence of the quadrature rule with increasing rif. The outcome is also shown 
in Fig. 3 as function of 

n tot = n t n Mcn it (82) 

and demonstrates an improvement by two orders of magnitude compared to the previous approach 
which fully evaluated the nine-dimensional integral by stochastic methods. Figure 4 shows a compar- 
ison with Gaussian integration which also gives fairly good results. 



This improvement now allows a much more precise determination of the coefficient (and, of 
course, also of the third order coefficient [37]). Table III contains the data for e^fi) from (3 = 4 
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Figure 4: (Color online) Comparison of deterministic integration routines for = 5) as a function 
of the number of integration points. The number of Monte Carlo calls {jimc) an d iterations (n%t) is 
kept fixed. An open symbol indicates a Monte Carlo result with inconsistent iterations. 



to /3 = 8 each with 12 iterations; the first 2 iterations were used for establishing the optimal grid 
while the following 10 were utilized for the statistics [denoted by n# = 12(2) in the following]. In 
addition to the classic VEGAS program (as in the previous test for 63) we also have used the VEGAS 
program from the CUBA library [38] which employs Sobol quasirandom numbers. This allowed to 
extend the range of inverse temperatures up to (3 = 10. Typical run times were about 1 day on 
a 2.4 GHz PC. It is seen that for all (3 there is agreement between the two data sets within the 
error bars. Despite larger statistics and higher accuracy in the deterministic integration the VEGAS 
(Cuba) routine returns larger errors which reflects our experience that the VEGAS (classic) error 
estimate often is too optimistic. This is also corroborated by the observation that at various /3-values 
the VEGAS (classic) results have an unacceptable large x 2 /Ndf indicating inconsistencies between 
different iterations within the given error bars. 
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Table III: Same as in Table I but for the fourth-order term e^{j3). The numerical results were obtained 
by a combination of deterministic and stochastic integration of the nine-dimensional integral (see text). 
Two different versions of the VEGAS program have been used: the classic one with pseudo-random 
numbers and the CUBA version with Sobol quasirandom numbers. The number of points in the 
deterministic tanh-sinh integration is denoted by n t . In the VEGAS (classic) evaluation nu = 12(2) 
iterations were used at each (3 value. Data marked by an asterisk have an unacceptable x 2 /Ndf 
(underlined) indicating that the iterations do not lead to a consistent error estimate. The last column 
gives the probability p that the error estimate for the VEGAS (Cuba) results is not reliable (p < 0.95 
is considered to be safe). 





VEGAS (classic): 


n M c = 4.7 x 10 5 


VEGAS (Cuba): 


n MC = 3 x 10 6 






n t = 23 




n t = 25 




-e 4 (/3) x 10 4 


X 2 /N DF 


-e 4 (/3) x 10 4 


P 


4.0 


0.4549 ( 6) 


0.637 


0.4563 (10) 


0.164 


4.5 


0.4828 ( 7) 


0.995 


0.4839 (11) 


0.157 


5.0 


0.5013 ( 8) * 


1.404 


0.5020 (12) 


0.170 


5.5 


0.5129 ( 8) 


1.087 


0.5136 (13) 


0.406 


6.0 


0.5193 ( 9) * 


1.739 


0.5209 (14) 


0.413 


6.5 


0.5239 ( 9) * 


1.488 


0.5254 (15) 


0.480 


7.0 


0.5271 (10) 


0.977 


0.5287 (16) 


0.534 


7.5 


0.5293 (10) * 


1.830 


0.5304 (18) 


0.588 


8.0 


0.5309 (11) * 


1.520 


0.5309 (17) 


0.646 


8.5 






0.5313 (19) 


0.387 


9.0 






0.5320 (19) 


0.355 


9.5 






0.5327 (20) 


0.483 


10.0 






0.5333 (19) 


0.553 
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But also for the VEGAS (Cuba) results the probability that the error is unreliable increases with 
the value of (3. This just reflects the fact that the cancellations inside the integrand are becoming 
more and more challenging at high (3. Fitting the VEGAS (Cuba) data with the asymptotic ansatz 
(64) yields 

e 4 = -0.5328(9) x 10~ 4 . (83) 
Data and best fit are shown in Fig. 5. The more general ansatz (65) leads to 

e 4 = -0.5330(7) x 10" 4 (84) 

with = 0.35(7). We therefore take 

e 4 = -0.533(1) x 10~ 4 (85) 

as our final result. 



0.54 H — 1 — 1 — 1 — 1 — 1 — 1 — L 




3456789 10 11 

Figure 5: (Color online) Same as in Fig. 2 but for the derivative of the fourth cumulant. The plotted 
data points are the VEGAS (Cuba) results from Table III. 
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C. A further step: e 5 



We have extended the BP approach to the calculation of the fifth-order coefficient es (/?) . Numerically 
this is much more challenging than the fourth-order calculation since these coefficients drop by roughly 
one order of magnitude in each order. This has to be achieved by cancellation in a 12-dimensional 
integral over a much more complicated integrand leading to much larger CPU times. 
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Figure 6: (Color online) Comparison of tanh-sinh and Gaussian integration for the derivative of the 
fifth cumulant at j3 = 4. Notation as in Fig. 4. 



Nevertheless the combination of deterministic integration and Monte Carlo integration leads to rea- 
sonable results. Figure 6 shows a slight advantage of the tanh-sinh integration rule compared to 
Gaussian integration. Of course, due to the more severe cancellations in the 12-dimensional integrand 
higher accuracy, i.e., a larger number of deterministic integration points is needed. At the same time 
the number of Monte Carlo calls cannot be as large as before to avoid excessive running times. 

Another numerical problem which already plagued the numerics for n = 4 (and to a much lesser 
extent n = 3) became more severe in the present case: due to round-off errors the Hadamard-Fisher 
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inequality (55) was not fulfilled exactly all the time: negative values down to 

xf£ = -3.1 x 1(T 9 (86) 

were recorded in double-precision arithmetic. Fortunately, this "digit-deficiency error" (see Ap- 
pendix B of Ref. [39]) does not affect the outcome of the Monte Carlo runs: checks have shown 
that e§ (/3 = 4) comes out the same whether the negative argument is set to zero or the absolute value 
of xhf is taken. In addition, the use of quadruple precision gives a consistent result (within error 
bars) but reduces the violation of the Hadamard-Fisher inequality considerably - at the price of a 
20-fold longer running time. 

Table IV: Same as in Table III but for the fifth-order term e^{0). For all deterministic numerical 
integrations fit = 25 integration points were used in the tanh-sinh integration routine. The Monte 
Carlo integrations were either done with the VEGAS (Cuba) program (timc = 1-5 x 10 5 ) or the 
classic VEGAS routine with umc = 7.9 x 10 4 ,n^ = 6(2) except for the data in boldface for which 
timc = 9.8 x W\n it = 5(2). 





VEGAS (classic) 




VEGAS (Cuba) 






-e 5 (/3) x 10 5 


X 2 /N DF 


-e 5 (P) x 10 5 


P 


4.0 


0.290 ( 4) 


0.240 


0.295 (10) 


0.369 


4.5 


0.337 ( 7) 


1.052 


0.317 (25) 


0.722 


5.0 


0.347 ( 6) 


0.537 


0.349 (18) 


0.353 


5.5 


0.365 (14) 


0.177 


0.330 (22) 


0.365 


6.0 


0.367 ( 7) 


0.287 


0.327 (26) 


0.657 


6.5 


0.361 ( 8) 


0.846 


0.370 (18)* 


0.956 


7.0 


0.365 (10) 


0.984 


0.394 (30) 


0.404 


7.5 


0.390 (13) 


1.296 


0.390 (42) 


0.329 




0.390 ( 9) 


0.592 






8.0 


0.366 (10) * 


2.514 


0.367 (35) 


0.326 




0.380 (15)* 


1.755 







The data are collected in Table IV and show that at high (3 it becomes more and more difficult to 
get consistent numerical results. Typical run times for each (3 value were about 1 month on a 3.0 GHz 
Xeon machine. With the Intel ifort compiler some loops could be vectorized leading to a reduction 
in CPU time by more than a factor of 2. If we exclude the data with x 2 /Ndf > 1-3 and p > 0.9 we 
obtain from a fit with zv 5 = 0.5 fixed 

e 5 = -0.378(4) x 10~ 5 . (87) 
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This is shown in Fig. 7 together with the corresponding values of 

(5) 3 

n tot = nt n Mcnu 



(88) 



for the different data from Table IV. It is not possible to determine the exponent v§ unambigously 
from the data which scatter too much. Taking a range of reasonable values for v§ we end up with 

e 5 = -0.38(2) x 1(T 5 (89) 

as final result for the fifth order energy coefficient. It is obvious that the given error is more an 
educated (and conservative) guess than a precise outcome of the fit. 




Figure 7: (Color online) Same as in Fig. 2 but for the derivative of the fifth cumulant. Data points 
with open triangles are from statistically inconsistent Monte Carlo iterations (see Table IV) and are 
not used in the fit. 
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VI. Conclusion and outlook 



We have shown that the Bogoliubov-Plechkov (BP) approach to calculate perturbative coefficients 
without diagrams works for the polaron problem (a simple field theory of electrons and phonons) 
if it is combined with several simple but crucial "tricks" to enhance the numerical feasibility and 
convergence. There is no indication that higher cumulants are "unbounded from below" as was 
reported in Ref. [40] in a much simpler anharmonic oscillator model [41]. It is worthwhile to point 
out the advantages and disadvantages of the BP approach compared to the standard perturbative 
method. 

While in the diagramatic approach a factorial increasing number of individual (zero-temperature) 
diagrams adds up to the final result, much fewer terms (moments) (see , e.g., Eqs. (27, 28)) must 
cancel inside the finite-temperature integral in the BP approach to obtain a result which is linear in (5 
so that the perturbative ground-state energy of the polaron can be determined. Of course, diagrams 
can be calculated exactly at zero temperature whereas in the BP approach the extrapolation (3 — ► oo 
must be performed numerically. We have demonstrated that by evaluating the derivative of the 
various cumulants, an exponential convergence to the zero-temperature limit can be exploited. Two 
new perturbative coefficients and es for the ground-state energy of a polaron have been obtained 
in this way and compared to results from Feynman's approximate treatment. 

It should be emphasized that the BP approach says nothing about the convergence of the perturbative 
series as it works in a fixed order. For the polaron case it is known that the ground-state energy is 
an analytic function of the coupling constant [42] but this is not necessary and systems where the 
perturbative expansion is known (or suspected) not to converge could be treated as well. Indeed, 
there is some hope that the methods which in the present work have been applied successfully for a 
simple nonrelativistic field theory may also be suited for relativistic field theories such as QED and 
QCD if these are formulated in the worldline formalism. Renormalization of the occuring divergencies 
is the main new challenge which is under investigation. 



Acknowledgement: Many thanks to Michael Spira who supplied his version of the classic VEGAS 
program and to Valery Markushin for help with compiler optimization which led to a considerable 
speed up of the calculations. I am also indebted to Dr. Plechko who informed me about his previous 
work in Ref. [8] and made some valuable remarks. 
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Appendix A: Energy coefficients from a quadratic trial action 

Here we briefly describe the results obtained with Feynman's variational method and with the best 
quadratic approximation [43]. Employing Jensen's inequality and working out the various path inte- 
gral averages one finds that the true ground-state energy is below the variational energy 



e < E t = n + v, 



where 



n = — 

2vr 



dE 



\nA(E) + 



1 



A(E) 



1 



V = - — 

/5F 



00 exp(— a) 
da 



o [/^V)] 1 / 2 



Here A(E) is the "profile function" which is related to the retardation function by 

sm 2 (Ea/2) 



poo 

A(E) = 1 + 8 / daf(a) 
Jo 

and n 2 (a) the "pseudotime" [44] given by 
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E 2 
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In Feynman's original work the retardation function is parametrized as 

f F {a) = Ce~ w ° 



(Al) 
(A2) 

(A3) 

(A4) 
(A5) 



which has the advantage that profile function, pseudotime and the kinetic term can be calculated 
analytically: 



A F (E) 
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(A6) 



Here v = ^/w 2 + 4C/w is used as parameter instead of the original strength C. Setting a = s 2 we 
thus have to minimize 
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(A7) 



where [45] 
b n (v) 



J_/-1/2 N 



ds e~ 



l-e~ 



t(» ](-«' 



2fe-l 
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(A8) 



For the actual calculation it is more convenient to introduce c = v /w — 1 = AC /w so that 



Ef(c,v) = -v 1 



1 



« - 



n=0 



(A9) 
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and to expand the parameters as 



c = c\a + C20? + . . . (A10) 
v = vq + v\a + V201 2 + . . . . (All) 

Including terms up to second order in a one finds vq = 3, c\ = 4/27, and Ep — ► —a — a 2 /81 — . . .. In 
higher orders the minimization always leads to linear equations for the coefficients c n ,v n so that they 
can be solved easily. With the help of a symbolic algebra system (such as MAPLE) the higher-order 
coefficients can then be evaluated in a straightforward manner and are given in Eqs. (14, 15). 

It should be noted that in lowest order also the retardation parameter w = 3 + 0(a) instead of 
w — > 1 as one would have expected naively. This is due to the wrong small-cr behavior in the ansatz 
(A5) for Feynman's retardation function and would be corrected by an "improved parametrization" 
[10] 

r , \ ol exi)(—wi<j) /A1 „x 

It is easy to check that both xi,wj — > l + 0(a) for small a. However, one can do even better by 
letting the functional form of the retardation function free. In this "best quadratic approximation" 
[43] one finds 

a exp(-cr) 

/best(CT) - 7-7= / M3/9 ( Aic V 

for which Eq. (Al2) is a convenient approximation since one knows that generally 

I^Vl^W^. (Al4) 
Indeed, inserting Uq(g) into the virial expression for the polaron ground state energy [43] 

a f°° exp(-cj) /3 \ , A1 _. 

one obtains -Evirial — ► —a for a — > 0, i.e., ek est = — 1. 

In second order we need the first-order change of the profile function and pseudotime 

A hest (E) = l + aa 1 {E) + a 2 a 2 {E) + ... (A16) 
^Lt(^) = ° + oiUxia) + a 2 U 2 {a) + . . . . (A17) 

From the connection (A3) between profile function and retardation function one finds 

4 f 00 , exp(-cr) sin 2 Ea/2 , , , 

and therefore from Eq. (A4) 

,exp(-<r') [°° sm 2 Eo/2sm 2 Eo' /2 
dE 



E 4 
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where ct< = min(cr, a'). It is possible to express the last integral exactly in terms of error functions 
and exponentials. However, for the calculation of the second-order energy it is better to plug this 
expression directly into the virial energy (A15) and expand [A t best( cr )]~ 1 ^ 2 U P ^° ^ rs ^ or der. 
Substituting a = s 2 ,cr' = s' 2 we then obtain 



2 ° 2 r i, e - r fiM_zfii e -.» + o ( Q ») (A20 , 



iTibest _ ^ 



10 

where s < = min(s,s'). Introducing polar coordinates s = r cos 4>,s' = rsincj) the integral with 
7r/4 < <fi < tt/2 can be combined with the one in which < (/) < ir/A and one obtains 

e !f st = -— / dr r 3 (3 - r 2 ) e~ r / # tan 2 (3 cos 2 <p - sin 2 M =-(—-—) . (A21) 

97T Jo JO V ' V 12 97T / 

Higher-order terms may be calculated numerically by using a delay-type equation for the pseudo- 
time which was found in the variational approximation for worldline QED and dubbed "variational 
Abraham-Lorentz equation" (VALE) [46]. It can be easily checked that the corresponding equation 
for the three-dimensional polaron case is 

2 , . d 2 u 2 (a) 4 f°° . 5V , w ,. 2a f°° , . exp(-cr') v . Ik . 

where 

A>, a') : = / u 2 cs t(^) " 2^bcs> + O - -^Lt (k " ^1) ( A2 3) 

is the delayed pseudotime (due to the phonon degrees of freedom which have been integrated out). 
Equation (A22) may be integrated with the boundary conditions /i 2 (0) = 0,/i 2 (0) = 1 to give 



^cst(-) = ° + ^ r da> {a - a>) r da " i e * p i x ( a '> ^ • (A24) 
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X bcst ( 

This gives an iterative scheme to calculate the perturbative terms (A17) for the pseudotime and 
eliminates the corresponding expansion (A16) for the profile function completely. Expanding in powers 
of a we obtain 

l_ r^r- ^ r,_///ex P (-a") 

3^7T 

Defining the delayed pseudotime of order n as 



U n (a) = ^ £ da' (a - a') J™ da" ' ^[jP Y n (a> , a') , n > 1 . (A25) 



X n (a', a") : = U n (a') - \ U n (a' + a") - X -U n (\a' - a"\) , n = 0, 1, . . . , (A26) 

the functions Y n are given by (for simplicity all arguments are suppressed) 

Y X Y~X 3XoUl Y-X 3 Aq^ + Xx^ 15X U 2 
Yi - A , Y2 - Ai - - — — , y 3 - A 2 - — + — — — , (A27) 



2 cr" 8 a' ; 

„ 3 X U 3 + X 1 U 2 + X 2 U 1 15 2X U 1 U 2 + X 1 U? 35 X Uf 
YA = X3 ~2 V' + ¥ ^ 16"^"- (A28) 



Once the perturbative terms U n (a) are known it is straightforward to calculate the energy coefficients 
g best^ n > j from the virial energy (A15) 
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with (again suppressing the argument a) 
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(A30) 
(A31) 
(A32) 



We have evaluated Eqs. (A25) - (A32) by numerical integration. This is a nontrivial task because 
of the square-root singularities at a = and the nonanalytic behavior of ^ 2 (\a — a'\). The first 



problem was solved by transforming to a = s , a = s , etc. , the second one by using the trapezoidal 
integration rule so that s = s' is precisely hit (and not integrated over). In addition, for the first 
three intervals of each integral a Newton-Cotes formula of open type [Eq. 25.4.21 in Ref. [47]] was 
employed in order to avoid evaluation of the various integrands at s = 0. While this cures the 
integrable singularities at the origin, it makes the treatment of the delay more problematic: in general 



U n (\a±a'\ 



\s 2 ±s' 2 \ 



is not in the tabulated values of U n (a = s 2 ) so that a three-term interpolation 



had to be used. In addition, the values of U n> ± for small a were determined from the (a = 0) limit 
of Eqs. (A22), (A23) 

o, n 2a f°° , exp(-o-) /A . 
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(A34) 



3-^/7? Jo V " 

with the same functions e n (<r) as used for calculating the energy coefficients. 

Although the trapezoidal (as well as the Newton-Cotes) integration rule is not very precise [it 
exhibits errors of 0(h 3 ) where h = s mSLX /N is the increment] it offers an additional advantage: the 
tabulation of U n (a = s 2 ) could be done step by step avoiding the time-consuming calculation of the 
integral over a' in Eq. (A24) for each value of a. Taking a max = 20 so that the retardation factor 
exp(— a) is sufficiently small at the upper limit of integration, we have achieved stable numerical results 
with N = 1000 — 1500. The numerical value of the second-order coefficent (A21) was confirmed with 
high accuracy (seven digits). 



Appendix B: Analytical results for the cumulants Ai and A2 

Here we calculate the cumulants \ n for n = 1, 2. In the first case an = t — t' = a and we have for the 
first moment 

rP r f . r°° exp(-(t-t')) 



a [p , r , . f°° , exp - 
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The remaining a integration is easily done by substituting s = a 2 . This gives 
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where erf(x) is the error function [47]. Thus we indeed have e± = — 1 for the first-order coefficient of 
the expansion of the ground-state energy in powers of the coupling constant. It is also seen that the 
subleading term in Ai is a constant which disappears if one calcualates the derivative of the cumulant 
with respect to 0: 
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The analytical calculation is more involved for n = 2. We start from Eq. (72) for the derivative of 
the second cumulant and substitute S = Si — £2 for the integration variable £i ;2 with £2,1 = 0— &2,i/2 
fixed. This gives 



d\ 2 

80 



4a 2 



vr jo 



ai da ^-^z^ r a dsfj au 



(B4) 



where s = (a\ + 02)/2 > r = (a\ — a 2) /I > 0. The explicit form (49) of a\ 2 may now be used to write 
the last integral in Eq. (B4) as 
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where the two parts correspond to the constant and linear behavior of 012, respectively, on the (S > 0) 
side of Fig. 1. We thus obtain 
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One sees that for — > 00 to (°i > ^2 , 0) since the relative times are bounded by the exponential 



retardation factors. In other words 
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and the corrections are of order exp(— 0). Putting 02 = t 2 a\, the o\ integration can be performed in 
the first term and an integration by parts in the second term gives 
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Finally a combination of partial integrations [to get rid of the arcsin(i) ] and integrals which MAPLE 
can do, leads to 

^ — ► a 2 [4 In (l + y/2) - 31n2 - V2~\ + O (e^ 3 ) . (BIO) 
We thus obtain the second-order coefficient of the ground-state energy as given in Eq. (7). 
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In our approach it is very important to know the precise way how Eq. (BIO) approaches the 
asymptotic value calculated above. The easiest way to find out is to differentiate Eq. (B6) again with 
respect to (3: 
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Consider first the contribution Iy. since to((3, 02, (3) = is the same as the upper limit of the 

integral, the latter vanishes so that 
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The substitution a 2 = (3s 2 gives 
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and in the limit (3 — > 00 the exponential factor forces s — > in all other terms [48]. Therefore we may 
expand these in powers of s, integrate term by term, and obtain 
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a 2 exp(— (3) 



For the contribution I 2 we use dto/d/3 = — Q(a\ + a 2 — (S)j ^Jo\o 2 so that 
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Using the variables of Eq. (48) we obtain 
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since the Jacobian of the transformation is 2. The substitutions s = (3(1 + u)/2, r = ssiiKp give 
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where sin^o('u) = (1 — + u ) • Again, for (3 — ► 00 the low-it behavior of the nonexponential part 

of the integrand determines the asymptotic behavior. We have 
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The <p integrals are elementary (see, e.g., Ref. [29], pp. 103, 104) and one obtains 
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Therefore 

/3^oo a 2 exp(-ff) 

/2 ( } 



is subasymptotic and after integration of Eq. (B14) with respect to (large) (3 we obtain 

/m 1 exp(-/3) 

e 2 (/?) — e,-^^-. (B21) 

Comparison with Eq. (B3) shows that this is the same functional approach to the asymptotic value 
as for the case n = 1 ; only the numerical coefficient is different. 



Appendix C: Tanh-sinh integration 

Here we briefly outline the "tanh-sinh integration" procedure proposed by Takahashi and Mori [33] 
and used in most of our deterministic calculations. For a one-dimensional integral over the interval 
x £ [— 1, +1] it is based on the transformation 

x = g(t) = tanh(Ksinht) i£ [— oo,+oo] (CI) 

g'(t) = 1 . — Kcosht (C2) 

cosh (Ksmhi) 

which has the effect that the transformed integrand g'{t) f(g(t)) vanishes at the boundaries along with 
all derivatives [for sufficiently well-behaved f(x)]. Therefore the Euler-Maclaurin summation formula 
[see, e.g. Ref. [47], Eq. 25.4.7] with a stepsize h does not get any (power) contributions from the 
endpoints and we have 



r+l p+oo k=+oo 

/ dxf(x)= / dtg'{t)f(g(t))^h ]T w k f(x k ) (C3) 

J - 1 k=-oo 



with 



x k = g(kh) = tanh [n smh(kh)] (C4) 

1 

cosh^ [k sinh(/c/i)] 



Wk = g'(kh) = ^Y~ r r~777vi K c osh(/c/i) . (C5) 



For large \k\ and fixed h we find 

x k — ► 1 — 2 exp (- K e |A;|fe ) , (C6) 
w k — > 2k exp (- Ke l fc l fe + \k\h) (C7) 

showing the "double-exponential" character of this transformation. 

Although the value k = n/2 has been reported to be optimal [35] we have found little difference 
in efficiency by taking 

k = 1 , (C8) 
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which is our choice in this work. In practice, the infinite sum in Eq. (C3) is finite since the weights 
Wk decrease rapidly with \k\ as seen in Eq. (C7). We use 

h\k\ < hk max = 3.4 (C9) 

as a cutoff so that x±k ma x = ±(1 — 2.01 x 10~ 13 ) and w±k ma x = 6.02 x 10~ 12 . The number of function 
calls then is 

n t = 2k max + 1 . (CIO) 
Conversely, if rit is chosen (as we do to estimate the run time in advance) the increment is given by 

fc= ™ (Cll) 
n t - 1 

It is straightforward to extend Eq. (C3) to an arbitrary integral as shown in Eq. (80) in the main 
text. 
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